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Abstract. The nucleus of the box-shaped galaxy NGC4013 
has been observed with the IRAM interferometer in the J=l-0 
and J=2— 1 lines of 12 CO. Our maps show the existence of a 
fast-rotating (130 kms -1 ) molecular gas disk of radius r~2" 
(HOpc). Several arguments support the existence of a bar po- 
tential in NGC4013. The figure-of-eight pattern of the major 
axis p-v plot, the ring-like distribution of gas, and the existence 
of gas emission at non-circular velocities are best accounted by 
a bar. 

We have also detected gas at high z distances from the plane 
(z^200-300pc=4"-5")- The latter component is related to a 
system of 4 Ha filaments of diffuse ionized gas that come out 
from the nucleus. The galactic fountain model seems the best to 
account for the Ha and CO filaments. Although the peanut dis- 
tortion can be spontaneously formed by a stellar bar in the disk, 
gas at high z might have been ejected after a nuclear starburst. 
The Ha filaments start in the plane of the disk at r^200pc(4"), 
and reach several Kpc height at r^600pc(10"), coinciding with 
the maximum peanut distortion where the strength of the restor- 
ing forces of the plane have a minimum. We have critically ex- 
amined other alternatives judged less probable: the existence of 
a CO warp (connected to the HI warping disk), the accretion of 
gas along stable inclined orbits and, finally, a vertical gas re- 
sponse near the resonances of the peanut (the latter is tested by 
numerical simulations). 

Although a link between the bar and the box-shaped bulge 
in NGC4013 is suggested we find noticeable differences be- 
tween the results of previous numerical simulations and the 
present observations. The discrepancy concerns the parame- 
ters of the bar generating the peanut. We see in NGC4013 the 
existence of a strong ILR region. The inclusion of a dissipa- 
tive component, which remains to be thoroughly studied, may 
change the evolution of the stellar peanut: although in simu- 
lations the peanut appears initially near a marginal ILR, the 
inflow of gas driven by the bar, can make two ILRs appear. 



Send offprint requests to: S.Garcia-Burillo 



Key words: galaxies: individual: NGC4013-galaxies: ISM- 
galaxies: kinematics and dynamics- ga 
laxies: spiral. 

1. Introduction 

N-body numerical simulations have established a link between 
the existence of boxy-shaped galaxy bulges and the presence 
of a bar instability in the disk (Combes and Sanders, 1981; 
Combes et al 1990, hereafter called C90; Raha et al 1991; 
Pfenniger and Friedli, 1991; Friedli and Martinet, 1993; Friedli 
and Benz, 1993). The bar can thicken, driven by vertical reso- 
nances, giving the bulge a boxy-shaped or peanut appearance, 
provided that the galaxy is seen at a high inclination (i>70°), 
and depending on the orientation of the major axis of the bar 
with respect to the line of nodes (measured by the azimuth <&): 
an edge-on bar ($=90°) adopts a boxy shape and the peanut 
feature appears for a large range of intermediate azimuths. The 
thickening is supported by star orbits that leave the galaxy 
plane and it occurs only in a privileged region along the bar 
where there is coincidence between the horizontal resonance 
ILR (0 6 

=0-ac/2) and its vertical counterpart (ilb=£l-i> z /2). Other me- 
chanisms, not related to the presence of a bar potential in the 
disk, can explain the formation of peanut bulges. Several au- 
thors support the merger/accretion interpretation (Binney and 
Petrou 1985, Shaw 1993, Shaw et al 1993). The occurrence of 
dissipationless collapse of an isolated system (Lima-Neto and 
Combes 1995), possibly in the presence of a bound companion 
(May et al (1985)) have also accounted for the existence of 
boxy shapes in galaxy bulges. 

In view of the large variety of scenarios it is necessary to 
envisage a detailed comparison between the different models 
and the observations for specific objects. There are very few 
examples where the connection between bars and box-peanut 
bulges have been investigated from an observational point of 
view. Bettoni and Galletta 1994 found photometric evidence of 
a bar in the peanut spiral NGC4442; Kuijken and Merrifield 
1995 studied a sample of edge-on spirals using Ha and optical 
lines as tracers of gas kinematics. However extinction effects 
can be severe in highly inclined disks where the column densi- 
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ties of gas and the associated dust can reach very high values. 
The use of macroscopically optically thin tracers such as CO 
can help us to get a more clear picture of the gas kinematics, 
especially in the central regions of these spirals. In particular, in 
this work we investigate the link between bars and box-peanut 
bulges using CO as a fair tracer of the molecular gas kinematics 
in the nuclear disk of a candidate galaxy: NGC 4013. 

The Sbc edge-on spiral NGC 40 13 has been identified as an 
extreme box-peanut bulge object, as seen in the optical pictures 
of van der Kruit and Searle (1982) and Jarvis (1986). Rand 
(1996) (hereafter referred as R96) have taken deep Ha im- 
ages of NGC 40 1 3 searching for extraplanar diffuse ionized gas 
(DIG), discovering an impressive set of filaments coming out 
the plane in the nuclear region. He suggests this is the signature 
of a nuclear outflow, after a starburst episode related with the 
distorted morphology of the nucleus. NGC 40 13 was mapped 
in the J=2-l and 1-0 lines of 12 CO with the IRAM 30m tele- 
scope (HPBW 13" and 21", respectively) by Gomez de Castro 
and Garcfa-Burillo, 1997 (hereafter referred as GCGB97). The 
CO disk extends to r=100" and it consists of a ring-like source 
(of radius r=30-40") and an unresolved fast-rotating nuclear 
disk. The high-velocity CO component has no HI counterpart 
(see Bottema, 1987, 1995 and 1996) and it is best explained 
by GCGB97 invoking the presence of a non-axisymmetric po- 
tential. Still the spatial resolution of the observations was in- 
sufficient to undertake a detailed analysis of gas kinematics in 
the nucleus of the galaxy. Furthermore the single-dish data sug- 
gested also the existence of out-of-the-plane gas that might be 
connected either with the vertical structure of the box-peanut 
bulge stable orbits, or the occurrence of an outflow and local 
ejections of material, as the DIG measurements of R96 seem 
to indicate. High-resolution (~3") interferometer CO data of 
NGC 40 13 are presented in this work, intending to study the 
distribution and kinematics of molecular gas in the nucleus of 
this galaxy and hopefully clarify the nature of the link between 
bars and box/peanut bulges. 

The study of dynamics of the gas in a peanut potential will 
deserve detailed numerical simulations involving both the stars 
and the gas, to be fully presented in a forthcoming publication 
(paper II). In this paper we limited ourselves to analyse whether 
there is a spontaneous and stable vertical response of the gas to 
a peanut potential, using a first run of numerical simulations. 

2. Observations 

The 12 CO(1-0) observations were made with the IRAM inter- 
ferometer of Plateau de Bure (see description in Guilloteau et 
al 1992), between December 1995 and March 1996. During 
less than one third of the observing period we observed simul- 
taneously the 2-1 line of 12 CO. We used the compact (CD) 4- 
antenna configurations, consisting of 18 baselines ranging from 
24 to 180 m. The declination of the source (S ^44°) allowed 
good UV coverage and made possible to synthesize an ellip- 
soidal beam of FWHM=3.4"x3.3" in the 1-0 line, adopting 
natural weighting and no taper on the visibilities (FWHM=1.2" 
xl" in the 2-1 line). The 12 CO(1-0) antenna half-power 



primary beam was 43" (22" in the 12 CO(2-l)). The pri- 
mary beam field was centered at a(J2000)=l l h 58 m 31.7 s , 
<5(J2000)=43°56' 

48.1". Assuming a distance of 11.6 Mpc (Bottema 1995), the 
43" primary beam corresponds to 2.4 kpc (l"=56pc). 

A total bandwidth of 430 MHz (1 120 and 560 kms" 1 for 
the 1-0 and 2-1 lines respectively), centered at v=839kms~ 1 
(LSR) and largely covering the range of velocities observed 
in NGC 4013 (±250kms-\ according to GCGB97), was ob- 
served with a resolution of 2.5 and 1.25 MHz (6.5 kms -1 and 
3.3 kms -1 at 1 15 and 230 GHz, respectively). The central 140 
MHz of this band was also observed at 1 15 GHz with a resolu- 
tion of 0.63 MHz (1.6 kins" 1 ). 

Data calibration was made in the standard way using the 
CLIC software package (Lucas 1992). The correlator was cal- 
ibrated every 20 min with a noise diode, and the RF passband 
once at the beginning of each observing run on 3C273. The 
relative phase of the antennas was checked every 20 min on 
the nearby quasars 1308+256 and 1156+295. The rms atmo- 
spheric phase fluctuations were typically between 10° and 25° 
at 115GHz. 

The data were then cleaned to yield channel maps. The rms 
noise in a 2.5 MHz-wide channel is 2.3mJy (ATs=0.018K)per 
3.35" 2 beam at 115GHz (the 1 -a noise is of 20mJy per beam 
in the 2-1 line or equivalently 0.35K). 

3. The CO maps 

In the following, we will denote by x and z the kinematical 
major and minor axis of NGC 4013 (x increasing to the NE 
and z to the NW). All coordinates will be referred to the dy- 
namical center (xc,zc) which has been determined as follows: 
we first fit the position angle of the galaxy (PA) together with 
the Zc position of the disk, applying a standard least-squares 
method to the integrated intensity map of fig2a. The inferred 
value for the orientation of the galaxy plane is PA = 66°±5°, 
in agreement with the value derived from optical imaging. Fi- 
nally both xc and the systemic velocity v sys are derived impos- 
ing the maximum symmetry on the the inner 10" region of ma- 
jor axis p-v plot taken at z=zc- We calculate (xc,zc)=(0,0)= 
a(J2000.0) = ll /l 58 m 31 s .36, 5(J2000.0) = 43°56'50".9. 
This position coincides within 1" with the radio-continuum 
sources detected at 6cm and 21cm (note that Bottema (1995) 
reports a wrong position for the radio-continuum source) as 
well as with the optical center determined by Palumbo et al 
1988. Similarly, the velocities (v) will be relative to v sys =&4-0± 
lOkms- 1 (LSR). 

Fig.l shows the 12 CO(1-0) velocity-channel maps ob- 
served with the interferometer (oriented along x and z axes). 
Emission appears from v=143 to -170 kms -1 , concentrated in 
a rotating edge-on disk whose vertical structure is marginally 
resolved with our 3.4" beam. However we notice 'out of 
the plane' gas excursions mostly in the bottom-left quad- 
rant (within the range Ax=(-5",25") and Az=(-3",-8")), vis- 
ible at certain velocities (within the range Av=(-20kms -1 ,- 
lOSkms- 1 )). 
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Fig. 1. 12 CO(1-0) velocity-channel maps observed with the Plateau de IRAM interferometer with a resolution (HPBW) of 
3.4"x3.1". Coordinates are with respect of the dynamical center (indicated by the cross). Velocity-channels range from 
v=150kms- 1 to v=-183kms- 1 by steps of e.Skms- 1 . Contour levels are -12, 12, 24, 40, 60, 80, 100 to 200mJybeam- 1 by 
steps of 30mJybeam _1 . 



Fig. 2. a(top): 12 CO(l-0) integrated intensity contours ob- 
served with the IRAM interferometer towards the center 
of NGC 401 3. x and z are offsets (in arcsec) with re- 
spect to the dynamical center C. Contours are 1 and 2 to 
24Jy.kms _1 beam _1 , by steps of 2Jykms~ 1 beam~ 1 . A zoomed 
view of figure a is displayed below, b(bottom): same as 
a) but for the 2-1 line of 12 CO. Contours are 4 to 15 by 
Uykms _1 beam _1 . 



At the derived dynamical center offset (0,0), CO emission 
displays a large spread of velocities: ±130kms _1 . As shown 
in Fig. 2a, which represents the 12 CO(1-0) velocity-integrated 
temperature map, the observed high-velocity gas at the cen- 
ter is linked with the presence of a conspicuous CO nuclear 
disk in the inner 100pc(2"). The nuclear disk shows an east- 
west asymmetry: the maximum of emission is located at x^2", 
i.e., eastwards with respect to the dynamical center. The re- 
markable asymmetry of the inner lOOpc of NGC 40 13 is more 
clearly shown by the higher resolution 2-1 data: Fig. 2b shows 
the velocity-integrated map in the 2-1 line which allows to re- 
solve the nuclear disk radially. The latter displays an asymmet- 
rical ring-like structure of radius r~1.7" (95 pc). 

The overall distribution of molecular gas in the disk shows 
a similar E-W asymmetry with respect to (0,0): CO is stronger 
and it is more extended in the eastern side, up to x=35", than 
in the western side of the disk, up to x=-25" (the single-dish 
data of GCGB97 already showed this asymmetry holds for 
radii 40" <r<100"). The same asymmetry is echoed by other 
star formation tracers (Ha and non-thermal radiocontinuum). 
A similar asymmetry in the CO distribution is present in the 
inner nucleus of our Galaxy (Bally et al 1987, 1988) and other 
spirals (see the case of NGC 891: Garcfa-Burillo and Guelin, 
1995). 

3.1. The mass of molecular gas 

If we take a CO to H 2 conversion factor of X=N(H 2 )/Ico=2.3 x 
10 20 cm- 2 K- 1 km- 1 s (Strong et al 1988), the H 2 mass derived 
from the 12 CO(1-0) interferometer map is M(H 2 )=3 x 10 8 Mq. 
For this we assumed the distance to be D=11.6 Mpc and in- 
tegrated the CO flux within a rectangular area of dimensions 
70"xl8", centered on the position (0,0). Including the mass 
of Helium, the total molecular gas mass in the Bure field is 
M 9as =M(H 2 +He)=4x 10 8 Mq. The nuclear disk mass content 
is derived to be close to M gas ^0.6 x 10 8 Mq (this is an upper 
limit as we integrated Ico within r=2" for all velocities, in- 



Fig.3. a(top panel) We overlay the Ha extraplanar emis- 
sion image of R96 showing the existence of 4 DIG fila- 
ments (grey scale) with the map of CO emission integrated 
within the velocity range Av=(0kms _1 ,-80kms _1 )(contours). 
We show an overlay of the 12 CO(1-0) velocity-channel 
maps at v=— 26kms~ 1 (b(top)),v=— 52kms~ 1 (c(middle)) and 
v=— 104kms _1 (d(bottom))(contours) with a K band image of 
the nucleus of NGC4013 taken from Shaw 1993 (gray scale 
contours). A sketch of the 4 extraplanar DIG Ha filaments is 
depicted in figures a-d (dashed lines). 



eluding the contribution of gas seen in projection close to (0,0) 
but located far from the nucleus in the plane of the galaxy). 

Taking into account that the shortest spacing measured 
by the interferometer is ^20m, we expect to filter scales 
<L>-20-25" at 115GHz (equivalently, <L>~10-15" at 
230GHz). We have derived the fraction of the single-dish 
12 CO(1-0) and (2-1) fluxes included in the Bure maps. The 
zero-spacing flux filtered out by the interferometer is kept very 
low in the 1-0 line: we detect nearly ^100% of the 30m-flux 
within the interferometer primary beam. On the contrary, the 
filtering is severe at 230GHz: we estimate that ~l/2 of the 30m 
flux is missing in the 2-1 primary beam. 

4. The vertical distribution of molecular gas 

Assuming the vertical distribution to be gaussian-like we have 
derived the thickness of the CO disk (Az) deconvolving the 
measured FWHM on the 12 CO(1-0) velocity-integrated inten- 
sity map of Fig. 2a by our 3.4" x 3.3" synthesized beam. The 
inferred Az shows no systematic radial trend along the ma- 
jor axis and it varies between 1.5" and 2.5", i.e. translated into 
spatial scales Az=80-130pc. These values are comparable with 
the thickness of the thin molecular gas disk in the Galaxy: 60- 
lOOpc (Bronfman et al 1988). 

The existence of CO emission at high z is clearly visible 
in the bottom left quadrant of Fig. 2a (see also Fig. 3a). Note 
that the lowest contour corresponds to 5x<ta, where a value 
of (7,4=0. 2Jybeam~ 1 kms~ 1 , has been obtained through the 
expression a A =j A adv, taking er=0.002Jybeam~ 1 and 
Av=100 

kms -1 . The mass of molecular gas at high z is 
M(H 2 )-1.5xl0 7 

Mq, (~5-10% of the total emission). The latter is a conser- 
vative lower limit as the CO-to-H 2 conversion factor might 
be significantly higher. The existence of vertical extensions 
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of molecular gas in the nucleus were already suggested by 
GCGB97 using coarse resolution single-dish data. 

The connection between the DIG filaments discovered by 
R96 and the above reported CO extraplanar emission is best il- 
lustrated at certain velocities. Figs.3b-c-d represent an overlay 
of a K-band image of NGC 4013, showing the peanut bulge 
distortion of the nucleus in the region (Ax, Az)= (±15", 
±10") with the CO-velocity channel maps at v=-26,-52 and 
-104kms _1 (Figs. 3b, 3c and 3d, respectively). A sketch of the 
DIG filaments of R96 is included for comparison. R96 inter- 
preted the impressive set of 4 DIG filaments above and below 
the nucleus of NGC 40 13 as the signature of a massive nuclear 
starburst. Supernova explosions might inject a huge amount of 
energy to the gas finally lifted to high z. At close sight, filament 
A is associated with CO emission breaking out of the plane 
(v=-26kms~ 1 ), although extraplanar gas appears to be more 
spread above the plane than the DIG feature (see v=-52kms~ 1 
and Fig 3a). Filament B is also associated with extraplanar CO 
(v=-104kms- 1 ). 

We have estimated the gravitational potential energy of the 
extraplanar molecular gas associated with filament A (the most 
impressive), using the inferred positions of CO emission above 
the mid-plane. Taking the two isothermal stellar components 
model of Jacobi and Kegel (1994), fitted on the luminosity dis- 
tribution out of the plane in NGC 40 13, the resultant gravita- 
tional energy per unit mass $ is 

$ = oflncosh(z/h) 

°t = (PlOl + P2cf )/ (pi + pi) 

where p, and <jj stand for the stellar mid-plane densities and 
velocity dispersions of the i=l,2 components and h is the fit- 
ted scale height (h=1.72kpc and a t =5 8. 3kms _1 for NGC 4013, 
according to Jacobi and Kegel (1994)). The average height of 
the CO extraplanar emission observed within the interval v=- 
13,-78kms _1 is z^270pc, which gives a total potential energy 
of 50x 10 51 ergs for M gas =4x 1O 6 M0. We emphasize that the 
above equation allows to calculate the total potential energy 
of the material at height z, and therefore should be taken as a 
lower limit to the total input energy that must be given to the 
molecular gas in the plane to reach the observed values: first we 
have not considered the kinetic energy of these gaseous struc- 
tures at present and we made no assumptions on the efficiency 
of the mechanism which is responsible for injecting the en- 
ergy to the gas in the plane. Whatever the nature of this energy 
source is, at least 50 type II supernovae are required (assuming 
that a type II supernova releases ~10 51 ergs to the ISM). This 
lower limit might transform into an order of magnitude higher 
number of required SN, once the radiative and dissipative en- 
ergy losses of the process are taken into account, henceforth 
suggesting that the nucleus of NGC 4013 might be experienc- 
ing a starburst episode. On the other hand, the restoring force 
of the galaxy plane is diminished near the position of the max- 
imum distortion of the potential and this limit can be slightly 
lowered. 



The FIR indices in NGC 40 13 do not indicate a massive 
starburst at the global scale of the disk: neither the normalised 
infrared fluxes, log(FIR/a^ aJ =-13.8 (where a ga i is the galaxy 
major axis in arcminutes) and log(FIR/Ls)=0.3, nor the IRAS 
colors (S25/Si2=1.5; Sioo/S25=30) are indicative of typical 
starbursts (see Huang et al, 1996 for discussion of a large 
galaxy sample). However there is evidence that supports a star- 
burst event in the nucleus: first, a huge fraction (~30%) of the 
total radio-continuum flux at 21cm is emitted by the nuclear 
disk; this fraction is the highest among the sample of galaxies 
published by Hummel et al 1991. Most noticeably, there is an 
association between the H-shape of the DIG filaments and ex- 
traplanar CO emission with the box-peanut appearance of the 
bulge. A similar set of filaments has been discovered by Rand 
in a galaxy also classified as a box-peanut: NGC 3079. Connec- 
tion between the nuclear starburst, the distorted stellar potential 
and the CO kinematics is discussed below (sections 5 and 6). 

HI observations (Bottema 1995) showed the existence of 
a highly warped disk in NGC4013. CO emission at high z 
present in our observations could come from regions of the 
outer warped disk seen in projection near the nucleus. However 
a comparison of velocity channel maps of atomic and molecu- 
lar gas leads us to reject this explanation: whereas CO emission 
at high z appears at negative velocities in the bottom left quad- 
rant (Figs. 2-3), the HI warp shows up in the top left quadrant 
(NE in Fig. 2 of Bottema 1995) within the same velocity range. 

5. Kinematics of molecular gas 

Fig. 4a shows the 12 CO(1-0) position-velocity diagram taken 
along the major axis. The high-inclination of the galaxy al- 
lows us to get the whole range of radial velocities in the plane. 
Molecular gas emission is spread in a romboid-like region 
where we distinguish three velocity components: 

- [I] The most outstanding component displays a large spread 
of velocities (v=±130kms~ 1 ) towards the dynamical center 
(referred as the straight ridge C) and it has no HI counter- 
part. The high-velocity gas corresponds to the central CO 
nuclear disk identified in Figs.2a-b and it extends from x=- 
3" to x=3". It is fully resolved in the 2-1 line where it ex- 
tends from x=-2" to 2"; note that dv/dr reaches a maxi- 
mum value of ~1000kms _1 kpc _1 in the inner 100 pc(2"). 

- [II] Part of the emission is concentrated in a straight ridge 
(R) which slowly drifts in velocity when we move along the 
major axis: it extends from (x=30", v=-80kms _1 ) to (x=- 
30", v=100kms _1 ). HI emission is detected in this region 
of the p-v plot according to Bottema 1987, 1995 and 1996. 
CO emission fills unevenly the p-v space between C and 
R. 

- [Ill] We detect gas at non-circular velocities (or velocities 
forbidden by circular rotation) in the inner ±5". This region 
(F) extends symmetrically on both sides of the nucleus:at 
x~-3.5", we measure CO emission up to v^-60kms _1 , 
whereas circular rotation would impose v>0 for x<0. The 
same applies for the offset x~3.5" where v~60kms _1 . 
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Fig. 4. a(left):The 12 CO(1-0) position-velocity diagram along the kinematical major axis, oriented along PA=66°. Gray contours 
correspond to 0.010, 0.020, 0.035 to O^OOJybeam" 1 by steps of 0.02Jybeam" 1 . A zoomed view on the inner r=5" region, where 
there is a CO nuclear disk showing high-velocity gas, is displayed on the right, b(right): same as a) but now for the 2-1 line of 
12 CO with a ~1" resolution in gray scale contours, linearly scaled from 0.035 to 0.200Jybeam _1 by steps of 0.025Jybeam _1 . 
Dashed contours stand for the same 12 CO(1-0) position-velocity diagram shown in a). Notation of kinematical components (C,R 
and F) are explained in text. 



Fig. 5. The radial variation of main frequencies of NGC 40 13 
(O, VL-kI2, Q+k/2 and £!+k/4) is plotted, assuming the epicyclic 
approximation and based on the CO rotation curve. The loci 
of the principal resonances are determined for a value of 
£!p=65kms~ 1 kpc~ 1 . 

5.7. The major axis p-v profile 

We have derived the CO-based rotation curve (v^) from the 
terminal velocities method, applied to the p-v major axis dia- 
gram. We have used the 2-1 data in the inner x=±2", the 1- 
data from 2" <x<25", and finally added the 30m data of 
GCGB97 for the outer disk (25" <x<100"). In the deriva- 
tion, we have taken into account the effect of channel smear- 
ing (Av=6.5kms~ 1 ) and assumed a typical cloud-cloud veloc- 
ity dispersion of ^lOkrns -1 . Finally we have assumed that v is 
translatable into v^, as circular- motions should be the main 
contributor to v (Sofue 1996; Garcfa-Burillo 1997). 

Note that the Hi-based rotation curve (v^,{) cannot account 
for the CO high-velocity gas component C. The scarcity of 
atomic hydrogen in the nucleus, together with the low reso- 
lution of HI observations can explain the differences between 
v^J t and v^\ We therefore conclude that the real rotation 
curve is certainly much steeper in the inner 20" than inferred 
using HI data: v^^HOkms -1 at r=110pc, this implies a dy- 
namical mass of M dyn — vxv^ ot / G— 2.8x10^Mq, assuming a 
spheroidal component in the nuclear region (a factor 0.6 lower 
in the case of a flat disk distribution). Therefore, we estimate 
M gas /Mdyn to be of 22% inside the nuclear disk. The latter 
ratio goes down monotonically to reach 10% at r=500pc. Al- 
though the dynamics is still dominated by the stars, the contri- 
bution of molecular gas to the total mass content of the nuclear 
disk is significantly high and the effects of gas self-gravitation 
might not be negligible. 

As stated above, the F component indicates the existence of 
emission at non-circular velocities. The high-resolution of the 
present observations precludes any explanation of F in terms of 
beam dilution: at x=±3.5" we are off by more than one synthe- 
sized beam from the center. 

5.2. A bar in NGC 401 3? 

The presence of the F component is more readily explained by 
invoking a deviation from axisymmetry in the inner (r^200pc) 
mass distribution of NGC 401 3. More precisely, we have ob- 



Fig. 6. The 12 CO(2-l) isovelocity contours (solid and dashed 
lines) in the nuclear disk of NGC 4013, linearly spaced from 
-lOOkms -1 to lOOkrns -1 by steps of 20kms _1 , overlaid on the 
corresponding integrated intensity map (fig.2b). 

servational evidence that the gas flow in this galaxy might be 
driven by a barred potential: 

- The major axis p-v diagram shown in Fig.4a displays the 
figure-of-eight shape typical of gas flowing along a bar 
(Binney et al 1991; Garcfa-Burillo and Guelin 1995; Kui- 
jken and Merrifield 1995). Components C and R, produc- 
ing a double-peaked line-of-sight velocity plot, would cor- 
respond to molecular gas populating two different families 
of orbits: x 2 inner orbits for C, and higher energy Xi orbits 
for R. The F component is best explained as the projection 
of the outer Xi orbits close to the cusped orbit. x 2 orbits 
develop between two inner Linblad resonances ILRs (the 
outer (oILR) and the inner (ilLR)). Xi orbits exist between 
the oILR and the corotation of the bar (COR). This particu- 
lar morphology of the p-v plot cannot be reproduced if we 
impose axisymmetry in the potential (Kuijken and Merri- 
field 1995). 

- Moreover we expect to see the imprint of a bar in the radial 
distribution of molecular gas as a reflect of secular evolu- 
tion. Gravitational torques induce a radial redistribution of 
gas which accumulates in rings between the ILRs, the 4: 1 
or UHR resonance and the outer Linblad resonance (OLR). 
As a result, the corotation region is progressively emptied 
of gas. The radial distribution of gas in NGC 4013 shows 
a compact central source of radius r^2" (C), a region rel- 
atively emptied of gas between C and R (r^ 10-30") and 
a secondary maximum towards r^50" (according to the 
single-dish data of GCGB97). This distribution, showing 
the existence of several rings, is well accounted by the bar 
hypothesis. 

- Fig. 6 shows the 12 CO(2-l) isovelocity contours in the nu- 
clear disk region. The z-distribution of molecular gas is spa- 
tially resolved at 1" resolution: we detect the presence of a 
slight velocity gradient along the minor axis of the galaxy 
and a westwards tilt of ^60° in the isovelocities of the nu- 
clear disk, suggesting that we are seeing in projection a 
non-axisymmetric gas distribution. 

Although the arguments enumerated above cannot be taken 
as a proof of the existence of a barred potential, it all points 
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out to this scenario as the simplest explanation. To explore on 
more firm grounds the consequences of this hypothesis, in the 
following section we try to locate the main resonances of the 
bar in the disk. 

5.3. Inferring limits on the bar pattern speed 

The standard first-order approach consists of deriving the prin- 
cipal frequencies of the disk from the fitted rotation curve 
(v^o?): Q, Q-k/2, Q+k/2 and Q+k/4 (see Fig.5). Assuming the 
epicyclic approximation, we can figure out the value of the bar 
pattern speed (Q p ) and consequently the position of the main 
resonances, based on observational and theoretical arguments. 
This procedure is not intended to provide an accurate fit of Q p , 
but to test if the basic results of numerical simulations are sup- 
ported by the present observations. 

As shown in Fig.5, the Q-k/2 curve presents a strong max- 
imum (~180kms _1 kpc _1 ) at r~3.5" and it goes monotoni- 
cally down to 1/10 of its peak value (^20kms -1 kpc -1 ) at 
r~20"; farther out, it stays quite flat until the edge of the opti- 
cal disk. Self-consistent numerical simulations of bars, includ- 
ing only the stellar component (C90), predict that the formation 
of the peanut occurs when the bar breaks up the plane near a 
marginal ILR (i.e., the bar pattern speed remains always close 
to the maximum of Q-k/2). If this is to apply in NGC 4013, we 
would require the bar to have an unlikely high pattern-speed 
(Q p >180 

kms _1 kpc _1 ) which also implies a marginal ILR at r^4". The 
latter is not compatible with the observations: the maximum 
distortion of the peanut bulge in NGC 4013 is at r~12". More- 
over, as we ignore a priori the orientation of the bar major axis 
along the line of sight, r~12" should be taken as a lower limit 
for the ILR. 

Most of the available numerical simulations have treated 
the appearance of the peanut distortion and its evolution in- 
cluding only the stars. Although the influence of a dissipa- 
tive component on the fate of the peanut remains to be studied 
thoroughly, a plausible evolutionary sequence can be advanced 
here. The peanut instability first sets in at a marginal ILR of 
the bar, but this should be taken as a starting point. The subse- 
quent inwards flux of gas towards the original ILR, driven by 
the bar's gravitational torque, can change the Q-k/2 curve in the 
inner region. The curve can become progressively steeper and 
finally two inner Linblad resonances (outer (oILR) and inner 
(ilLR)) are bound to appear. The stellar bar itself is expected 
to slightly slow down in the process. The peanut distortion is 
the relic of an old marginal ILR; however a strong ILR region 
will develop in the course of time. Any reasonable value of Q p 
in NGC 40 13 implies we have two ILRs at present: an inspec- 
tion of fig.5 leads to that conclusion unavoidably. Therefore, 
the existence of a peanut instability and a strong ILR region in 
the nucleus may not be in contradiction, but should be taken 
as a result of secular evolution. Numerical simulations to be 
developed in paper II will allow us to test this scenario . 

Fitting the bar pattern speed is beyond the scope of the 
present paper, however a lower limit on Q p can be tentatively 



set. The radial distribution of molecular gas (with a relative de- 
pression or hole between r=10" and r=30") suggest that coro- 
tation cannot be at a radius larger than r^50", implying that 
Q p >60-70kms _1 kpc _1 . The latter implies the following loci 
of principal resonances: i iILR >2", i oILR <12", i C or <50", 
*uhr <70" and v lr <100". Corotation of the bar is kept 
well inside the optical disk. This is expected for bars in spi- 
rals of early or intermediate Hubble types such as NGC 4013, 
classified as Sbc (see Combes and Elmegreen 1993). Secondly, 
if Q p >60-70kms _1 kpc _1 , the oILR appears inside the old 
marginal ILR (at 12"), as proposed lines above. Finally the 
major axis distribution of CO showing a maximum at r~50" 
is well accounted for if a ring forms near the UHR resonance 
(Whr <70"). 

6. A link between the gas filaments and the bar? 

In a barred and peanut-shaped potential, the gas in principle 
tends to follow the periodic orbits, and is dragged through ver- 
tical resonances out of the plane, like the stars. However, the 
gas is dissipative, on a short time-scale; cloud-cloud collisions 
occur with a characteristic collision time of 10 7 yr. While the 
peanut and perpendicular elevation time scale is a long process, 
of time-scale almost 10 9 yr (C90). Therefore, the gas compo- 
nent looses its disordered kinetic energy, and in particular in 
the z-direction, rather quickly. The gas settles down in a very 
thin disk, irrespective of the peanut-shape potential. 

In order to study if there can be a stable and spontaneous 
vertical response of gas in a peanut potential, we have per- 
formed self-consistent numerical simulations involving gas and 
stars. In this first run, the gas mass fraction is kept very low 
on purpose (~2%). The gas clouds are considered as sticky 
particles and we neglected the effects of self-gravitation. No 
effect of star formation on the gas dynamics is considered ei- 
ther. The code used is FFT-based particle-mesh, with the gas 
treated as sticky particles (Combes & Gerin 1985). The useful 
grid is 128x128x64, and the suppression of the Fourier im- 
ages is done through the algorithm of James (1977). The run in- 
cluded 150k stellar particles and 40k gas particles. A bar forms 
spontaneously in the stellar component, since the bulge to disk 
mass ratio is 1/3, and there is no other spheroidal component 
that could stabilize the disk. After 2xl0 9 yr, the bar has de- 
velopped a characteristic peanut shape (see contours of Fig. 7). 
The main result is that, at any time, the gas disk remains very 
thin (~ 200pc at half-power), and is never perturbed by the 
peanut/shaped potential or only in a very transient manner (see 
Fig. 8). 

The vertical thickening of the bar and the subsequent for- 
mation of the peanut can be caused by instabilities associated 
with resonances between the bar motion and the vertical oscil- 
lations of the stars(C90; see also Binney(1981)); this effect has 
also been attributed to the buckling or fire-house instability by 
Raha et al 1991. The bar instability first acts to increase the 
eccentricity of stellar orbits and align their principal axes; this 
causes the buckling instability, precisely about the vertical res- 
onance region, which increases vertical velocity dispersion and 
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Fig. 7. Contours of stellar projected density of the bar seen 
edge-on, after 4 Gyr of simulation. The scale is in kpc. 

thickness. Once the bar is thickened, the nature of the stellar or- 
bits in the peanut can be described with orbits trapped around 
the 2:2: 1 periodic orbit family: these correspond to the vertical 
Lindblad resonance, where in the frame of the bar, the particles 
just perform two z-oscillations in one turn (SI — v z /2 = Oj>, 
the bar pattern speed). This happens to occur in the region of 
the in-plane inner Lindblad resonance (where — k/2 = 
and therefore the resonant orbits have the 3-D shape of a ba- 
nana (projected into an ellipse in the plane (C90, Pfenniger & 
Friedli 1991). A detailed orbital study has emphasized the bi- 
furcations of the main periodic orbits x\ onto the banana and 
anti-banana orbit families (Pfenniger & Friedli 1991). These 
periodic orbits, in a strong bar, very often possess loops, and 
the gas is not likely to follow them because of dissipation. This 
may explain the small propensity of the gas to be elevated ver- 
tically, at least for directly rotating orbits. For retrograde or- 
bits, on the contrary, the orbit family X4 bifurcates due to ver- 
tical instability into the anomalous orbits (corresponding this 
time to a 1/1 resonance, or one vertical oscillation for one turn, 
Pfenniger & Friedli 1991). This family does not possess loops, 
and the gas can be maintained more easily on these (Friedli 
& Benz 1993). These retrograde orbits could be populated in 
particular during an accretion event, where gas clouds arrive 
with a given angular momentum, un-correlated with that of the 
galaxy. The case of NGC 128 is a remarkable example. It is a 
peanut-shaped galaxy seen edge-on, where an inclined gaseous 
disk is observed to counter-rotate with respect to the stars (Em- 
sellem & Arsenault 1997). Although this galaxy has apparently 
no sign of perturbed morphology, it has a companion nearby 
(NGC 127) which could have provided the retrograde gas. 

These stable retrograde orbits are characteristic of a tum- 
bling triaxial potential. Already van Albada et al (1982) had re- 
marked that in the rotating frame of a triaxial object, the Cori- 
olis forces on retrograde particles produce a torque, which is 
opposite to and compensates the torque of the restoring forces 
towards the plane. The latter would have made the retrograde 
orbits to precess, and since they differentially precess with ra- 
dius, the clouds' collisions and dissipation would have forced 
the gas to settle in the main plane. The Coriolis forces, there- 
fore, stabilise the gas into the inclined retrograde orbits. These 
orbits are not particularly related to the peanut shape, but both 
phenomena indicate the presence of bars. 

In NGC 4013, the presence of a prodigious warp (Bottema 
et al 1987) suggests that a large amount of gas has recently 
been accreted with un-related angular momentum. Part of it 
could have gone towards the center. If spontaneously, directly- 
rotating gas will not stay at large height above the plane due 
to the peanut shape of the potential, it is possible that accreted 
gas with a different angular momentum, not aligned with the 
principal axis, could be trapped in banana orbits, or more likely 
in retrograde anomalous orbits. More simulations are needed, 



taking into account the accretion origin of the gas, and also its 
larger mass fraction and self-gravity (see paper II). 

An alternative mechanism to explain the presence of gas at 
high altitude above the plane, is the galactic fountain effect, due 
to massive star formation. This could explain the morphology 
of the Ha filaments, that appear much higher above the plane 
than the neutral gas, traced by the CO emission in the center. 
The fact that the filaments seem to coincide with the peanut 
region might not be a coincidence: the gas is easier to expel 
when the restoring force of the plane is less, i.e. in the thicker 
stellar plane. 

Extraplanar neutral+ionized gas and dust have been found 
in another edge-on spiral: NGC 891. The existence of a thick 
molecular disk, containing M(H 2 )^1O 8 M0 up to z^lkpc, was 
first established by Garcfa-Burillo et al (1992) using single dish 
CO observations. R96 also reports the detection of an extended 
DIG halo with some prominent Ha filaments in this galaxy. 
Recently Howk and Savage (1997) have discovered the dusty 
counterpart of the thick CO disk: the HST and WIYN opti- 
cal BVR images of NGC 891 show hundreds of dust filaments 
lying far from the mid-plane. The derived neutral gas mass 
content of these absorbing structures closely agrees with the 
CO-based estimates. Although some of the extraplanar dust 
features are interpreted as supernova-driven galactic chimneys, 
other are less clearly linked with highly energetic phenomena 
in the disk. In either cases, the existence of an extended neu- 
tral+ionized gas with a set of dust filaments is related to the 
high star formation efficiency in the disk (Garcfa-Burillo et al 
1992; R96; Howk and Savage 1997). 

7. Summary and Conclusions 

We have observed with the IRAM interferometer the emission 
of the 1-0 and 2-1 lines of 12 CO in the nucleus of NGC 4013, 
an edge-on Sbc spiral possessing a box-shaped bulge, with spa- 
tial resolutions of ^3.3" and ~1.2", respectively. 

Our maps show the presence of a distinct fast-rotating (v rot 
~130kms _1 ) nuclear disk of radius r^2" (lOOpc) and gas mass 
M gas ~0.6x 10 8 Mq. The high-velocity component (absent in 
the HI map) is accompanied by gas emission at non-circular 
velocities within ±3" from the dynamical center, indicating 
that the gas flow is driven by a non-axisymmetrical potential. 
An analysis of the gas kinematics, derived from the 12 CO(1-0) 
major axis position-velocity plot and the 12 CO(2-l) isoveloc- 
ities map, supports a bar model for NGC 4013. The observed 
ring-like distribution of molecular gas (at r^2"(100pc) and 
r~50" (2.7kpc); the outer ring inferred from the 30m map) is 
interpreted as the imprint of a bar. 

Although a link between the bar and the box-shaped bulge 
in NGC 4013 is suggested, there are apparent discrepancies be- 
tween the results of numerical simulations and the model pro- 
posed here. Disagreement concerns the basic parameters of the 
bar generating the peanut. In self-consistent simulations of the 
star component the peanut forms near the locus of a marginal 
ILR of the bar. Instead, the derived Q-k/2 curve in NGC 40 13 
leaves little doubt of the existence of a strong ILR response, 
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Fig. 8. Particle plots at the same epoch (4 Gyr) of the stars (left) and gas (right). Up is the face-on view, and bottom the edge-on 
view. In each frame, the maximum radius is 25kpc. 



irrespective of the chosen pattern speed of the bar (fl p ). The 
nuclear disk would be located in the vicinity of the ilLR, with 
no gaseous ring counterpart in the oILR. The present observa- 
tions suggest that the inclusion of a dissipative component in 
simulations might probably change the evolution of the stel- 
lar peanut: although the latter appears near a marginal ILR, the 
inflow of gas driven by the bar, makes two ILRs appear and ac- 
celerates the secular evolution of the gaseous disk in the inner 
~lkpc region of the galaxy. Numerical simulations, involving 
both the stars and the gas in a peanut potential, will analyse the 
response of the gas in the disk of the galaxy fully testing the 
bar hypothesis (see paper II). 

We analysed the vertical distribution of molecular gas, 
showing that, although the bulk of CO emission comes 
from a thin disk (with deconvolved thickness FWHM^80- 
130pc), there is a non-negligible amount of molecular gas 
(M(H 2 )=1.5x 10 7 Mq) at large z distances from the plane 
(z^200-300pc). The close relationship between the DIG fila- 
ments seen in Ha coming out of the plane and the presence 
of molecular gas emission, suggests that both share a common 
origin: gas ejected by a massive nuclear starburst. 

A preliminary run of simulations has restricted to study the 
vertical response of the gas to a peanut potential, that sponta- 
neously forms in a disk of stars. Gas clouds are treated as test 
particles and we neglect here the effects of star formation and 
self-gravity in the gas dynamics. Due to its dissipative nature, 
the gas forms a thin disk very quickly. Contrary to the stars, 
the gas cannot be maintained at high altitude above the galaxy 
plane along stable orbits. Cloud-cloud collisions make impos- 
sible the population of banana and antibabana self-intersecting 
orbits. Moreover, it remains to be studied if gas clouds can pop- 
ulate the vertical bifurcation of the retrograde X4 family in a 
peanut, after an accretion episode (see paper II). However this 
mechanism is unlikely to hold for NGC 4013, where the bulk 
of the gas in the disk is directly rotating. 

Among the different explanations for the gas at high z- 
inclined resonant orbits connected to the peanut, gas accretion 
in the course of an interaction and, finally, the galactic foun- 
tain model- the latter seems the best to account for the Ha 
and CO filaments. Although the peanut distortion formed in 
the stars comes from a bar in the disk (the presence of the latter 
being suggested by the observed CO kinematics), gas is being 
ejected in the nucleus after a bar driven starburst. The filaments 
come from the inner 200pc(4") and reach a height of several 
Kpc, coinciding with the maximum peanut distortion where the 
strength of the restoring forces of the plane is diminished. 
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